STABILITY INVESTIGATIONS OF AIRFOIL FLOW BY GLOBAL ANALYSIS 
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As the result of global, non-parallel flow stability analysis 
the single value of the disturbance growth-rate and respective 
frequency is obtained. This complex value characterizes the 
stability of the whole flow configuration and is not referred 
to any particular flow pattern. The global analysis assures 
that all the flow elements (wake, boundary and shear layer) 
are taken into account. The physical phenomena connected 
with the wake instability are properely reproduced by the 
global analysis. This enhance the investigations of instability 
of any 2-D flows, including ones in which the boundary layer 
instability effects are known to be of dominating importance. 
Assuming fully 2-D disturbance form, the global linear stabi- 
lity problem is formulated. The system of partial differential 
equations is solved for the eigenvalues and eigenvectors. The 
equations, written in the pure stream function formulation, 
are discretized via FDM using a curvilinear coordinate sy- 
stem. The complex eigenvalues and corresponding eigenvec- 
tors are evaluated by an iterative method. The investigations 
performed for various Reynolds numbers emphasise that the 
wake instability develops into the Karman vortex street. This 
phenomenon is shown to be connected with the first mode 
obtained from the non-parallel flow stability analysis. The 
higher modes are reflecting different physical phenomena as 
for example Tollmien-Schiichting waves, originating in the 
boundary layer and having the tendency to emerge as insta- 
bilities for the growing Reynolds number. The investigations 
are carried out for a circular cylinder, oblong ellipsis and air- 
foil. It is shown that the onset of the wake instability, the 
waves in the boundary layer, the shear layer instability are 
different solutions of the same eigenvalue problem, formula- 
ted using the non-parallel theory. The analysis offers large 
potential possibilities as the generalization of methods used 
till now for the stability analysis. 

Introduction 

The boundary layer stability analysis based on the so- 
lution of the Orr-Sommerfeld equation is a useful tool for 
practical analysis of the laminar- turbulent transition. The 
only competing method is based on purely empirical formu- 
las, characterized most often by the shape parameter. 

It is widely accepted that infinitely small disturbances, 
although amplified according to linear stability theory are 
not able to onset the laminar-turbulent transition unless the 
amplification reaches some value so a factor has to be intro- 
duced to correct the results of the analysis. The method 
has been develop to match the results of the empirical an 
theoretical investigations. 

The laminar-turbulent transition is usually preceded by 
the Tollmien-Schiichting waves. Several receptivity experi- 
ments (Morkovin [4]) were provided to understand the phe- 


nomena of the Tollmien-Schiichting waves generation. It is 
commonly accepted that Tollmien-Schiichting waves are ge- 
nerated by an external source of disturbance (as for example 
acoustic excitation) and that the non-parallel or non-uniform 
effects enhance the feedback between the wave and the ex- 
citation. These non-parallel and non-uniform effects are the 
viscous boundary layer growth, the change of the surface 
curvature and variation of the surface static pressure. The 
growth of the boundary layer is evident near the leading edge 
of the blunt body, change of the surface curvature causes 
the non-parallelity of the flow, surface static pressure chan- 
ges significantly in the separation region. It is characteristic 
that these three problems were studied separately. Goldstein 
[1] solved analytically the problem of evolution of Tollmien- 
Schiichting waves near the leading edge. The influence of 
sudden change of the geometry was investigated by Gold- 
stein [2] and Ruban [3] . 

These investigations have one common feature - the as- 
sumption of slow variation of the flow in the streamwise direc- 
tion as necessary condition for weakly non-parallel analysis. 

Elliptic nature of the Navier-Stokes equation describing 
the flow suggest that the phenomena in all these regions are 
not independent and influence each other. The question ari- 
ses if interactions of the leading edge geometry, boundary 
layer and wake can be described by a single theory. The 
natural choice is to drop the parallel flow assumption and 
to treat the flow in all these regions as a whole. The con- 
sequence is the attempt to use the non-parallel flow, global 
stability analysis. The non-parailel theory was succesfuly 
used to study the wake instability [10, 9, 8, 7, 12]. There 
are no theoretical limitations to apply this analysis also to 
various geometries, as for example the airfoil. Because the 
assumptions of the non-parallel theory is a generalization of 
the classical parallel flow analysis, one can expect that this 
method is adequate not only for determination of the wake 
instability. The instability of the boundary and shear layer 
must be reflected in the eigenvalue solutions of the problem. 

Governing equations 

Linear stability theory is concerned with the development 
in time and space of infinitesimal perturbations around a gi- 
ven basic flow. If this basic flow is assumed to be paral- 
lel, the classical theory of parallel shear flow stability can 
be applied. This method has been also successfully used for 
nearly parallel flows for which the multi pie- scale method, ad- 
opting the concept of ’’slow” variation of flow parameters in 
one direction, is valid. In general, non-parallel case only the 
two-dimensional theory taking into account the non-parallel 
effects is adequate. The equations of this theory are briefly 
presented here. 

The problem was solved in the pure (Lagrangian) stream 
function finite difference formulation. This formulation, not 


very common in the Navier- Stokes equations solvers, offers 
certain advantages for the eigenvalue analysis. The primitive 
variables formulation ([9]) results in much larger matrices. 
Although the eigenvalues are equal for velocities and pressure 
one has to dead with the full system. This difference in size 
is even more evident because the matrix entries are complex 
for the eigenvalue analysis. 

The unsteady incompressible Navier* Stokes equations writ- 
ten in the stream function formulation take the form: 

[| + (Vx^.v-1a]a^o (1) 

$ = iWs (2) 

We assume that the stream function 0(x t y,f) is a sum of a 
steady part V>(x, y) and the unsteady disturbance ^'(x, y,t) : 

=$(x,v) + i> , (x,y,t) (3) 

The disturbance value is assumed to be small compared to 

the stream function value. Introducing equation (3) into (1) 
we obtain the nonlinear equation: 

d 1 

|^+(Vx^)-V — — A] = 0 (4) 

Assuming a small disturbance allows the linearization of the 
equation (4) i.e. we ignore the terms containing W>*) 7 . In 
the disturbance equation we separate the time and space de- 
pendence: 

^'(x T y,f) = p(x,y)e~ iXt (5) 

where 

A = 7T (St + i<j) (6) 

Introducing the above relationship into (4) results in the li- 
near partial differential equation: 

iAA<£- (V x rfr) • VA(,p — (V x p) • VA^ + “-A 2 <^ = 0 (7) 

He 

The fundamental difference between this equation and the 
Orr-Sommerfeld one, which is derived in similar manner as- 
suming the disturbance form as: 


The symbol | denotes the covariant derivative of the function. 
For further specialized metric tensor coefficients 


Sn = a *(0 0 ({, 7 ) 

927 = 0 3 (v)9((,v) (12) 


only g(( , 77) and its first order derivatives g ^ and g tf) have to 
be calculated for any transformation. 

Reynolds number Re and Strouhal number St are expres- 
sed as: 




(13) 


Discretized, equation (11) can be written as: 


[A - XB)p = 0 (14) 

and represents the generalized eigenvalue problem. 

For the eigenvalue calculations complex numbers can be split 
into real and imaginary parts so that only the real arithmetic 
haw to be applied. Then the two parts of equation (2.10) may 
be written: 

Aipr “ A r Bp r + X^Bpi = 0 

Api - A iBp r + A r J Bpi = 0 (15) 


Solution 


Numerical discretization and meshes 

The discretization of the Navier-Stokes equations (11) 
and disturbance equation (12) is accomplished using the fi- 
nite difference method. In both cases the thirteen-point sten- 
cil was used. The accuracy of the derivatives for such a stencil 
is maximum 0 (h 2 ) for the fourth order terms. The unsteady 
version contains implicit stepping in time. 

For all the calculations the orthogonal O-type mesh ob- 
tained by the conformal mapping is applied. The Karman- 
Trefftz transformation was used for the airfoil calculations. 
The metric coefficients (13) are expressed analytically by 
means of symbolic manipulation program to assure the ma- 
ximum accuracy. 


= y)' i{ax - 0t) ( 8 ) 

is that, while Orr-Sommerfeld equation is an ordinary diffe- 
rential equation, equation (5) is a partial differential equa- 
tion. This means different methods of solution and numerical 
problems encountered for the two cases. 

To solve the problem for an arbitrary flow geometry the 
curvilinear body fitted coordinate system should be used for 
the solution of the equation (1) and (7). For orthogonal 
metric the following relations are valid: 

9x, = o , g' 1 = 0 , i ^ j (9) 

hence equations (1) and (7) can be written as: 

0 °) 


Boundary Conditions 

For the steady Navier-Stokes equation solution the follo- 
wing boundary conditions are used: 

y> = 0 , 0,n = 0 on the body (16) 


= 0 , V>,n =Vpot,n in the far field (17) 

The collocation of the vorticity transport equation is made 
only for the outflow. For the inflow the Dirichlet boundary 
condition with the value of the potential flow solution is ta- 
ken. The boundary conditions for the disturbance equation 
(12) are: 

on the body (18) 

in the farfield (19) 


P = 0 , p, n = 0 


Du) 

Ti 


= o 


Dp 

Ti 


= 0 



The Dirichlet boundary condition (zero disturbance) is in- 
troduced for the inflow. The introduction of the convective 
boundary conditions appears to be an important factor of 
improving the numerical accuracy, especially for the steady 
and unsteady flow calculations. 


Solution of the eigenvalue problem 

In any eigenvalue problem the question arises whether all 
the eigenvalues are sought or whether determination of only 
one or few is satisfactory. Solving similar problem Zebib and 
Kim et al. [10,11] applied the QZ type decomposition from 
the standard libraries. The advantage of finding all of the 
eigenvalues is that no guess values have to be made. For rela- 
tively small matrix size, resulting from the use of the spectral 
method or crude FDM meshes this procedure is acceptable 
and was used in our earlier investigations [7]. Jackson ap- 
plied for the un symmetrical, complex generalized eigenvalue 
problem, appearing in the non-parallel flow stability theory 
the inverse iteration method [9], This concept is also adopted 
in our present investigations. The eigenvalue, closest to the 
guess value and the related eigenvector are both determined 
at the same time. Till now it is the only realistic method for 
very large equation systems. 

The following equations explain the principle steps of this 
method. Applying the Newton-Raphson method to equation 
(14) we obtain 


(A - X WB)(v>W + V"’) - dX^B^ n) = 0 (20) 


which can be written as: 

(A - AWf?)r7 (n+1) = (21) 


where the normalization is performed as follows: 
^( n+1 ) = 

and 

(O, = S ir 


( 22 ) 

(23) 


denotes a unit vector. The correction of is calculated 
from: 


dA( n+J > = 


1 

{e r ) T T}( n + 1 ) 


(24) 


The iteration process involves the repeated solution of the 
equation (21), normalization of the eigenvector and correc- 
tion of the eigenvalue. This process continues until conver- 
gence of the eigenvector and eigenvalue is achieved. The pro- 
cedure, which consists of LU decomposition at each step with 
a quadratic rate of convergence, was replaced by a method 
using only one LU decomposition. The convergence is then 
only linear but the back-substitution time is significantly re- 
duced compared to the decomposition time, justifying many 
iteration steps: 


(A - A (25) 

The scheme is found to be convergent to the eigenvalue clo- 
sest to A 0 and to produce the appropriate eigenvector. 

.Numerical results 

The linear stability analysis consist of two steps. First 
the steady solution of the Navier- Stokes equations has to be 
found. In practice both, the steady and unsteady solution of 








Figure 1: Steady flow solutions for the circular cylinder flow 

the Navier-Stokes equations was performed. The unsteady 
one served as the reference data for the comparison to the re- 
sults of eigenvalue analysis. It is characteristic that obtaining 
of the unsteady solution near the critical Reynolds number 
is difficult. For symmetrical flow some external forcing has 
to be introduced. The response of the flow field is dependent 
on the way the disturbance is introduced. The nearly neutral 
stability of the flow caused that the influence of the distur- 
bance dominates the flow even after a long time. In this 
case the purely numerical aspects of the computation are of 
much greater significance. Also unsymmetrical flows near the 
critical Reynolds number requires a lot of CPU time to be- 
come fully unstable. The flow patterns of initial periods are 
different from the "fully developed” unsteady ones (Fig. 12). 
Near the critical Reynolds number such patterns can persist 
over a long time requiring significant amount of CPU time 
to obtain the real periodic state. Some codes fail to carry 
out the calculations long enough in time and due to unphysi- 
cal boundary conditions the solution breaks down when the 
vorticity reaches the outflow boundary. The unsteady simu- 
lation for the Reynolds number higher than the critical one 
is easier. For this reason always the higher Reynolds number 
unsteady solutions were taken for the comparison with the 
stability analysis. 

In the linear stability theory the Navier-Stokes equati- 
ons are linearized about a steady flow. The quality of the 
steady solution has then the direct influence on the eigen- 
value analysis. The accuracy of the solution is the best for 
the circular cylinder flow and is decreasing for the ellipsis 



Figure 2: The growth-rate and the Strouhal number for the 
circular cylinder flow. 


and airfoil flow where leading and trailing edge can cause 
numerical problems even for meshes generated by the con- 
formed mapping. In case of limited computer resources it is 
satisfactory for the numerical simulation of the flow to use 
relatively crude mesh Bpacing on central, upper and lower 
parts of the airfoil. In this case the gradients of the quan- 
tities along the boundary layer are not very large. For the 
eigenvalue analysis however, also the fine discretization in 
this direction is very important. The attempt to detect the 
Tollmien-Schlichting waves necessities at least several tenth 
of points for one period preserving also the fine discretiza- 
tion in the radial direction. The compromise for these two 
contradictory requirements was partly obtained by calcula- 
tion of the steady solution on one mesh and interpolation of 
the result on another mesh, more suitable for the stability 
calculations. 

The eigenval ue solution wy ^culate d^ foi^^he exter nal 
flow around the circular cylinder, ellipsis andean airfoil. The 
circular cylinder served as the source of reference data, for 
the validation of the program because a lot of numerical and 
experimental results is avaiable. The only existing results 
for non-parallel analysis are the circular cylinder results [10, 
9]. The flow around the ellipsis was investigated to analyze 
different eigenmodes. The modes characterized by higher 
frequency are clearly appearing for high Rey golds nu mber s . 
Because of the extremely long wake for Re > 200, causing 
several numerical difficulties such an analysis could not be 
carried out for the circular cylinder. Finally the N AC A 4412 
airfoil flow for a = 0° and a = 15° was shown to examine the 
potential possibilities connected with the eigenvalue analysis 
of this geometry. 

Circular cylinder results 

For the symmetrical flow around cylinders it is always, 








Figure 3: Real (a) and imaginary (b) part of the eigenvector. 


theoretically, possible to obtain a steady-state solution, even 
above the critical Rey nolds number. The streamlines pat- 
terns obtained for the steady flow around a circular cylinder 
are shown in Fig.l. These results served as the input data 
for the eigenvalue analysis. The guess value for the Strouhal 
number is 0.12 and the growth-rate 0. The result of the cal- 
culation consist of the complex eigenvalue for each Reynolds 
number together with a complex eigenvector. The growth- 
rate and the corresponding frequency as the function of the 
Reynolds number is shown in Fig. 2. Some results of our pre- 
vious investigations using the QZ method are also plotted. 
The results of these calculations are compared with those 
obtained by Zebib [10], which uses the non-parallel analysis 
in the spectral stream function formulation together with a 
full-matrix eigenvalue solver of a QZ-type. For the inverse 
iteration method, used in our computations, the critical va- 
lues are Rtc = 46.23 and St c = 0*1345. 

The read and imaginary part of an eigenvector for the_ in- 
creasing Reynolds number is depicted in Fig. 3. Over a wide 
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Figure 4: Eigenvector velocities (imaginary part): (a) below 
Rt c (Re = 40) (b) above Re c (Re = 50). 


range of Reynolds numbers the eigenvector (disturbance) pat- 
terns are very similar, showing the physical aspects of the 
phenomena to be already present in flows of fairly small 
Reynolds number. The increase in Reynolds number allows 
these modes to cross the zero-growth-rate line and emerge as 
instabilities. The problem arises if there is any difference in 
eigenvector patterns bellow and above the critical Reynolds 
number. It is known from the parallel flow stability analy- 
sis that the wake stability is governed by its characteristics 
in the vicinity of the rear stagnation point. Careful study 
of the eigenvector values near the cylinder shows (Fig.4) the 
difference in the disturbance patterns above the Re c . This 
enhance the onset of the Karman vortex street. 

To evaluate how realistic are the obtained eigenvalue so- 
lutions the disturbance is summed with the steady-state solu- 
tion for Re = 90. As the reference the unsteady flow simula- 
tion for Re=l00 is taken (Fig. 5). The same periodic patterns 
are present in both pictures. This proves that for the cylinder 
flow instability the non-linear effects are not significant. 

Ellipsis flow 

Following the approach for the circular cylinder flow the 
elliptic cylinder was analyzed. It is known from experiments 
and non-parallel flow stability analysis of Jackson, performed 
for the bodies with different cross-sections that the proper 
scaling of Strouhal number is based on the dimension per- 
pendicular to the main flow direction. For such a scaling its 
value is not much different for various shape of the cylinder. 
The critical Reynolds number reflects also the overall shape 
of the body. The relation between the axis ratio of the el- 
lipsis and the critical Reynolds number was studied earlier 
[8]. For the oblong ellipsis situated parallel to the flow di- 
rection the critical Reynolds number is increasing while the 
slope of the growth-rate curve becomes smaller, comparing 
to the circular cylinder results. As can be expected the Kar- 
man vortex street mode results differ only slightly from ones 
obtained for the circular cylinder. The eigenvector patterns, 
growth-rate and frequency relations for increasing Reynolds 
numbers are similar to the circular cylinder ones. The in- 
teresting results are obtained also for the Reynolds number 
higher than the critical one. We assume that the steady flow 
solution coincides with the real one in the boundary layer and 
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Figure 5: Karman vortex street (a) superposition of the di- 
sturbance and steady solution, Re = 90 (b) unsteady simu- 
lation, Re = 100 
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Figure 6: Higher mode eigenvector (real part) for the 1:5 
ellipsis flow. Re _= 200 — • 


the shear layer near the body, even for the Reynolds number 
higher than the critical one. The justification for such an 
assumption are the experimental investigations of Kourta et 
al. [13] and Unai and Rockwell [14] in the higher Reynolds 
number range the Karman vortices are formed not directly 
behind the cylinder. Between the cylinder and the vortex 
street a dead fluid zone is found, bounded by two nearly 
parallel shear layers. As the Reynolds number increases the 
length of the dead-fluid zone decreases and the location of 
the first instability waves in the shear layer moves upstream. 
According to the results of the parallel flow stability analysis 
the unsteady behavior of the fluid is governed by the flow in 
direct neighborhood of the body. This conclusion allows us 
to cut the steady solution and limit the computational do- 
main. The fact that the length of the wake, obtained as the 
steady-state solution of the Navier- Stokes equations exceeds 
the assumed "infinity” distance (the wake end is outside the 
computational domain) is in context of the eigenvalue ana- 
lysis not relevant. 

This steady flow solution is was used as the base for the 
eigenvalue analysis. The assumed guess frequency is higher 
than for the Karman vortex mode. The result of the higher 
mode analysis is depicted in Fig. 7 and 8. The growth-rate 
is a function of both Reynolds number and mode, so that 
different modes are preferentially amplified as the Reynolds 
number increases. In Fig. 7 the growth-rate and the Strouhal 
number for higher mode is depicted together with the first 
one for the ellipsis having the axis ratio 1:5. The temporal 
evolution of the waves is shown in Fig. 8. The amplitude of 
the wave is raising in the direction of the separation. The 
waves on the upper and lower surface of the ellipsis are shif- 
ted in phase as the result of superposition of the symmetric 
pattern of disturbances and antisymmetric stream function. 
The characteristic patterns for all higher modes investigated 
are the family of branches of disturbance streamlines having 
sequentially positive and negative values. Each branch is 
ended with a ceil located in the vicinity of the maximum 
velocity gradients in the boundary or shear layer. The eigen- 
vector patterns should be analyzed in connection with the 
steady flow solution. The disturbance is added to it to ob- 
tain the unsteady flow. In the steady solution two regions 
can be distinguished - "soft” part where the stream function 



Figure 7: The growth-rate and the Strouhal number for the 
for 1:5 ellipsis flow 

values are small, containing the boundary layer, separation 
and wake region and the "stiff" part where the stream func- 
tion values are large in comparison to the disturbance. It is 
obvious that when adding the disturbance and steady state 
solution only the "soft” part is "modulated” while the "stiff" 
one is practically not influenced (Fig.8). For thjs reason the 
considerations concerning the eigenvector patterns outside 
the "soft” region have very limited practical meaning. This 
conclusion is confirmed by numerical calculations, showing 
that the "soft" regions of the eigenvector are related to the 
growth-rate and frequency value. The rest of the field is more 
likely influenced by numerical aspects of the computations. 

For the Blasius profile instability the Tollmien-Schlichting 
wave length is approximately six times larger than the boun- 
dary layer thickness. Since the boundary layer on the ellipsis 
is relatively thick for the range of the Reynolds numbers ap- 
plied in the calculations the detected Tollmien-Schlichting 
waves are also long. The shorter ones, for higher Reynolds 
numbers require much finer meshes, especially in the circum- 
ferential direction. The eigenvector cells, located on the ellip- 
sis surface near the leading edge are shorter (in the circumfe- 
rential direction) than the ones in the separation region. For 
a given constant frequency which is the same for the whole 
field it can mean only that the wave propagates slower near 
the leading edge and faster in the separation region. The 
propagation along the shear layer of the wake has approxi- 
mately constant velocity. All the found eigenvalues for the 
Tollmien-Schlichting mode were damped ones. The question 
arises if the Tollmien-Schlichting wave, considered globally, 
in the boundary layer and propagating further along shear 
layer can become amplified without external excitation. The 
growth-rate is raising with the increasing Reynolds number 
and one can expect that the higher mode wave will become 
only slightly damped or even amplified for the high enough 
Reynolds number. 

For any flow around the cylinder exist many eigenmo- 
des. In practice near any given frequency exist an eigenvalue, 
mostly with such an low growth-rate that it is unlikely that 
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Figure 9: Steady flow solutions - NACA 4412, a = 15° 

it can emerge as the instability. Similar conclusions can be 
drawn on base of the Kim [11] results. 

Different eigenvectors can be classified into at least two 
groups. One characteristic eigenvector pattern is connected 
with the onset of the Karman vortex street. Fig.3 shows this 
mode for the circular cylinder. Jackson [9] has shown the 
same patterns. Similar mode was detected by Karniadakis et 
al. [5] who investigated the flow around the circular cylinder 
placed in the channel bounded with two parallel plates. This 
mode is called there the central mode and dominates for the 
cylinder placed near the symmetry axis. Moving the cylinder 
toward the wall causes switching to the "wall mode" which is 
related to the Tollmien-Schlichting waves. For the external 
flow around the cylinder the "wall” mode forms similar cells 
located however on the body and in the shear layer. 



Figure 8: Tollmien-Schlichting waves - temporal evolution 
for the 1:5 ellipsis flow, Re = 200 


_The airfoil flow 

The another cylinder flow which was considered is the 
airfoil flow. As the example geometry the NACA4412 airfoil 
is taken. Two different angles of attack were considered. 
For q = 15° the stall is evident and the regular Karman 
vortex street appears for high enough Reynolds number. The 
numerical simulation of such a flow was performed by Shutz 
[6]. For a = 0° dominating phenomena take place in the 
boundary and shear layer. 

First the steady flow solution has been found (Fig. 9). The 
character of the steady flow solution for a = 15° is different 
from the circular cylinder one. (Fig-1, Fig.9). While for the 
circular cylinder the wake consist ot two bubles, there is only 
one for the airfoil flow. 

The eigenvalue analysis gave the fastest growing mode 
(Fig.10). 

For a = 15° the flow becomes unstable at Re = 335. The 
eigenvector patterns are in this case also very similar to ones 
for the circular cylinder (Fig-11 ) . In Fig. 13 the comparison 











Figure 10: The growth-rate and the Strouhal number for the 
airfoil flow 

between the real part of the eigenvector for Re - 100 and 
Re ~ 500 is shown* The value of the disturbance is growing 
with the flow direction for both cases. It is normalized, so 
the disturbance reaches the same maximum, located in the 
vicinity of the outflow boundary. Because for Re = 100 
(Fig. 13) the growth-rate is negative the disturbance will be 
damped after a long enough time. The flow for Re = 500 
is unstable. The disturbance is growing both in time and in 
the flow direction. The characteristic feature for the higher 
Reynolds numbers flows is the much larger amplitudes of the 
disturbance in the wake close to the airfoil. 

To compare the obtained eigenvalue analysis results with 
the real flow patterns the unsteady simulation was used. The 
simulation was performed for Re = 1000. The early stages 
of unsteady simulation exhibit patterns significantly different 
from the "fully developed" ones (Fig. 12). This discrepancy 
is even greater in the neighborhood of the critical value. For 
this reason to compare with the eigenvalue analysis one pe- 
riod was taken after long enough time (t = 56.8 to t = 64.0). 
Earlier periods are "spoiled” by the initial flow development. 
The comparison of the flow patterns for Re = 600 (eigen- 
value analysis) and Re = 1000 (unsteady simulation) show 
very good qualitative agreement. All the mechanisms of the 
vortex shedding are properly reproduced. This fact is one 
more proof that the Karman vortex street, especially near 
the body has the linear character. 

For the angle of attack equal 0° till Re = 800 exists no 
separation on the airfoil. The higher mode solution forms two 
row of cells (Fig.15) which are close to the airfoil only near the 
leading edge. When added to the steady flow solution only 
the shear layer behind the airfoil is effected (Fig.16). The 
flow is stable because the growth-rate is negative, but if it 
becomes unstable it is the Kelvin-Helmholz type of instability 
of the shear layer. For increasing Reynolds numbers the cells 
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Figure 11: Real (a) and imaginary (b) part of the eigenvector 
- airfoil flow, a = 15° 


are moving closer to the airfoil. The disturbances form now 
cells attaching the airfoil and forming the "wall” mode. The 
boundary layer is now "modulated” in the way simlar to the 
ellipsis flow. For o = 0° the Karman vortex street mode also 
exists, although it is strongly damped for the small Reynolds 
numbers. 


Conclusions 

It was shown that non- parallel flow stability analysis is a 
method most suitable for determination of the wake flow in- 
stability. Several examples , calculated for different Reynolds 
numbers and geometries ranging from circular cylinder to the 
airfoil with the angle of attack, show that the method is a 
general tool for prediction of the wake instability. It is of 
advantage of this method, comparing to other numerical ap- 




Figure 15: Higher mode solution for the NACA 4412 airfoil 
a = 0°, a) Re = 300, b) Re = 900 



Figure 16: Superposition of the steady solution and higher 
mode disturbance for the NACA 4412 airfoil, Re — 300 

proaches, that the critical Reynolds numbers and respective 
frequencies are determined more precisely. The method is 
able to handle the unsymmetrical wake flow. In this calcula- 
tions the superiority of the stream function formulation and 
iterative determination of the eigenvalue has been proved. 

Using the same method higher modes were investigated 
for the ellipsis and airfoil flow. Although the investigations 
had a preliminary character it can be concluded, that the 
results obtained differ significantly from the first mode so- 
lution. The higher mode disturbance patterns, obtained for 
the ellipsis flow, added to the steady solution appear to be 
the Tollmien-Schlichting wave originating in the boundary 
layer. The wave propagates further along the mixing layer. 
For the airfoil flow these types of modes were found but also 
another instability phenomena axe present in the eigenvalue 
solutions. The investigation of large spectrum of eigenmodes 
is even more difficult because some instability phenomena 
are smoothly "switching” to another ones. The result of the 
higher-mode analysis gives the qualitative insight into the 
stability problem. The limitation of the method on present 
state of its development is not the formulation but the num- 
merical approach. These difficulties we hope to overcome in 
our future investigations. 


We belive that the method presented here will enable 
the stability analysis of any flow as a whole, without brea- 
king it into pieces or restricting considerations to single type 
and that all instability phenomena are reflected in the non- 
parallel flow eigenvalue solutions. 
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